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Can we do better than Hybrid Monte Carlo in Lattice QCD ? 

M.E. Bcrbenni-Bitsch a , A.P. Gottlob a , S. Meyer a * and M. Piitz a + 

a Fachbereich Physik - Theoretische Physik, Universitat Kaiserslautern, 
D-67663 Kaiserslautern, Germany 

The Hybrid Monte Carlo algorithm for the simulation of QCD with dynamical staggered fermions is compared 
with Kramers equation algorithm. We find substantially different autocorrelation times for local and nonlocal 
observables. The calculations have been performed on the parallel computer CRAY T3D. 



1. INTRODUCTION 

Overcoming critical slowing down in Monte 
Carlo simulations of lattice field theories is im- 
perative to approach the continuum limit. Col- 
lective mode algorithms have improved the qual- 
ity of numerical studies for systems with bosonic 
degrees of freedom substantially during the last 
years. Besides many other things the important 
role played by autocorrelation functions has been 
recognized and in particular strong finite size ef- 
fects present in the decay of autocorrelation func- 
tions have been observed. 

On the other hand it is fair to say, that the 
improvements of fermionic simulation algorithms 
for lattice field theories since the introduction of 
the Hybrid Monte Carlo (HMC) algorithm ten 
years ago have not been overwhelming. The the- 
oretical analysis of the HMC algorithm has been 
possible in the free field case while the practi- 
cal analysis to compare the performance of vari- 
ants of the HMC algorithm and other dynamical 
fermion algorithms for lattice QCD has only be 
done for Wilson fermions &0]. 

The work reported here is a first step to obtain 
reliable estimates for integrated autocorrelation 
times of different operators in lattice QCD with 
dynamical staggered fermions using two different 
simulation algorithms. 
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2. AUTOCORRELATION FUNCTIONS 

In lattice quantum field theory one evaluates 
expectation values of observables A(<f>) 
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A{cj>) 



with S{4>) the action and 4> x a field on a lattice A 
by a dynamic Monte Carlo algorithm. 

The dynamics of the numerical algorithm is a 
stochastic process which is realized on a computer 
by a Markov process. The transition probability 
matrix P{<j) — > 4>') leaves the equilibrium distri- 
bution invariant. 

Detailed balance is equivalent to the selfadjoint- 
ness of P as an operator on L 2 (y) with real spec- 
trum (X m in, ^max), and the stationarity condition 
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= 1 • e 



may be read as an eigenvalue equation with eigen- 
value A = 1 and eigenvector e~ s ^\ All other 
eigenvalues A have the upper bound 1. 

The exponential autocorrelation time T exp pa- 
rameterizes the gap between A = 1 and the sub- 
leading (unwanted) modes. The expected error 
a a of an observable A((f>) is 
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for n > T exp 



with the autocorrelation function 
C AA (t) = ((A(i)-A)(A(i + t)-A)) 
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and the integrated autocorrelation time 
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The performance of dynamic Monte Carlo al- 
gorithms for a finite system with linear size L and 
correlation length £ is described by empirical dy- 
namical scaling laws 

T mt ,A ~ min(L, £) Zi "t. A . 

3. GENERALIZED HMC ALGORITHM 

Introduce a set of "fictitious momenta" it and 
a Hamiltonian 
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then the HMC alternates two Markov steps: mo- 
mentum refreshment and Monte Carlo molecular 
dynamics which contains molecular dynamics and 
a global Metropolis step to correct for the dis- 
cretization errors in Hamilton's equations. 

The first step is momentum refreshment, where 
momenta tt are replaced by new values chosen at 
random from a Gaussian distribution. 

Horowitz Q suggested a generalized momen- 
tum refreshment (Kramers equation) by adding a 
Gaussian noise 



tt' = e- 1&T ■ tt + a/I - e- 2 ^ T • r\ 
with 

P(V) = |e-" 2 / 2 

where < 7 < 00. For 7 = 00 one obtains full 
momentum refreshment, while 7 = gives exact 
generalization of the second order Langevin algo- 
rithm. Detailed balance will be satisfied if the 
new phase space configuration is accepted with 
probability 

P[(0,tt) -> (0', tt')] -min(l,e H(0 ' 7r) - ffW '^' ) ). 

4. NUMERICAL RESULTS 

All our results are for lattice QCD with four 
flavours of dynamical staggered fermions and 
gauge group SU(2). 



In the HMC algorithm there are two free pa- 
rameters: the integration step size St and the 
trajectory length tq = n <5r, where n is the num- 
ber of integration steps. St must be adjusted 
to keep Pace reasonably large. Optimal choice 
for To is less clear. In the Gaussian model with 
L>^ the dynamical critical exponent z is z = 2 
for To = constant, and z = 1 if To oc £. For 
constant To the dynamical critical exponent z be- 
comes z — 2 for £ L. 

For dynamical fermions with m = 0.1 we mea- 
sured autocorrelation times for the plaquette, 
Polyakov loop, and the chiral order parameter 
with trajectory length 1/4 < To < 2 on various 
lattices with sample sizes of ~ 1000 Tint- 
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Figure 1. Autocorrelation functions of the pla- 
quette data. 



In Fig. 1 we show the autocorrelation function 
of the plaquette on a lattice of size 8 4 for both 
algorithms. For the Kramers equation algorithm 
k = 12 (see below) and on the horizontal axis t is 
given in units of molecular dynamics time. The 
straight lines are fixed exponentials to guide the 
eye. 

For the Kramers equation algorithms the sec- 
ond exponential in the decay of the autocorrela- 
tion function seemed to be stronger coupled. This 
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observation has been found for local as well as 
nonlocal observables. 

4.1. Generalized HMC algorithm 

Through the generalized momentum refresh- 
ment a new tunable parameter 7 is introduced. 

Detailed balance requires that the momenta 
must have their sign flipped after every rejected 
step. 

The Monte Carlo molecular dynamics step can 
be performed k times. For dynamical fermions 
with m = 0.1 we measure autocorrelation times 
Tint, A for the plaquette, Polyakov loop and the 
chiral order parameter with 0.03 < St < 0.12, 
and k = 4,8 and 12. 
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Figure 2. Comparison of the HMC and Kramers. 



For more details we refer to [pj . 

In Fig. 2 we compare the cost in molecular 
dynamics time for the Polyakov loop operator as 
a function of stepsize for both algorithms. Data 
with k = 4 and k = 8 are included and compared 
with HMC using a trajectory length tq = 1. We 
consider the pattern of this plot similar to the 
findings of Q in their figure 2. 

4.2. HMC on the CRAY T3D 

The computational demands of our study let us 
move to the massively parallel processing CRAY 
T3D with a peak rate of 150 Mflops in 64 bit 



arithmetic per node. We decided to implement 
the shared memory routines because they show 
less overhead than Parallel Virtuell Machine or 
Message Passing Interface. We found the ex- 
pected speedup within a few percents for 4 PEs 
up to 128 PEs on a 16 4 lattice. For our present 
implementation of the HMC algorithm 5 PEs on 
a T3D are equivalent to one YMP processor indi- 
cating the improvment of the performance to cost 
ratio @. 

5. CONCLUSIONS 

We find that reasonable high statistics are ab- 
solutely necessary to obtain reliable error bars for 
the integrated autocorrelation times. We consider 
the different decay behaviour of the integrated au- 
tocorrelation functions for the observables stud- 
ied so far as an indication for different efficiency. 
Up to now we did not find a competitive set of 
parameters for the Kramers equation algorithm. 
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